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ABSTRACT 


The scope of this thesis is to study the stability of two ships in close proximity 
towing. Unlike previous studies in the past, the lateral dynamics of both ships are 
included in the formulation. The equations of motion of the system consist of the sway 
and yaw motions of the two ships and a control law for the leading ship. An eigenvalue 
stability analysis of the coupled system confirms the results that are obtained through 
numerical simulations. It is shown that it is possible for the system to be unstable even 
though the classical criteria for towing stability are satisfied. A series of parametric 
studies is conducted in order to analyze the sensitivity of the system for different towline 
lengths, tension, and control time constant. 
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I. INTRODUCTION 


The stability of ship towing is the primary concern in both naval architecture and 
maneuvering control. The main problem consists of two individual tasks, first building a 
mathematical model of the two towed ships and second analyzing the directional stability 
of the combined system. 

Traditionally, ship towing in the open ocean has been performed under relatively 
large separation distances between the two vessels; i.e., large towlines. Only in the case 
of barge towing in rivers and canals, has close proximity towing been applied. Over the 
last few years, however, interest on ocean close proximity towing has been resurfaced. 
The Office of Naval Research is interested in close proximity towing for its Advanced 
Hull forms Program, in particular small waterplane area twin hull (SWATH) ships and 
their variations (such as the SLICE hull). Studies of the SEA LANCE project at the Naval 
Postgraduate School (Total Ship Systems Engineering Program) have also shown several 
benefits associated with close proximity towing. One of the main benefits of close 
proximity towing in military applications is the ability to separate a combatant ship from 
a main part of its payload. 

Towing is not without its problems, however. Directional stability under tow and 
seakeeping are two main areas of concern. This thesis will concentrate on the problem of 
directional stability. Previous studies on directional stability of ship towing were 
performed by Abkowitz in 1964 who developed the characteristic equation for single 
body towing, and by Bemitsas and others who developed the criteria for stability of 
Abkowitz’s 4 //! order characteristic equation. In these studies the stability criteria was 
based solely on the trailing ship, which means that the dynamics of the leading ship were 
neglected. In fact, the leading ship was approximated by a point mass moving on a 
straight line under constant forward speed. Our concern is that although this may be valid 
for long towlines, it may not be a valid approximation for close proximity towing, and the 
existing stability criteria may be inadequate which is the reason and the justification of 
this thesis. 
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In this study as a leading ship we will use the SLICE vessel, which is 105 ft, and 
180 tons, and as a trailing ship the SWATH ship Kaimalino which is 89 ft, and 217 tons. 

The Slice is a high speed variant of the Swath technology. It has 4 underwater 
hulls instead of two. Attached to each hull is a strut that extends up to support the main 
body. Figure 1 shows the profde view of the Slice. [2] 



The SSP Kaimalino was the world’s first high performance open ocean Swath 
ship. It consists of two parallel torpedo-like hulls. Attached to the hulls are two 
streamlined struts. The struts extend above the water surface and support the main body. 
The Kaimalino also has stabilizing fins attached near the aft end of each hull. Figure 2 
shows a profile of the SSP Kaimalino. [2] 



Figure 2. The profile of the SSP KAIMAFINO [From: [2]] 


2 




































The general approach adopted in this thesis is as follows: 

1. First we will use the maneuvering equations of motion in the horizontal plane for 
both leading and trailing ship. 

2. Coupling between the two ships will be provided through the towline. 
Hydrodynamic coupling arising from radiation and diffraction effects will be 
neglected. This can be easily incorporated into the analysis once the effect of such 
hydrodynamic coupling on the hydrodynamic coefficients is established. 

3. The response of the coupled system will be analyzed by both simulation and an 
eigenvalue analysis. It is hoped that this analysis will provide regions of 
directional stability, and can therefore lead to a design methodology for rational 
towing system parameter selection based on their stability properties. 

The hydrodynamic coefficients of both ships used in this study are provided from 
reference [2], It should be emphasized that our procedure will be general enough so that it 
can accommodate different or additional hydrodynamic coefficients. 
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II. PROBLEM FORMULATION 


A. EQUATIONS OF MOTION 

The surge, sway and yaw maneuvering equations of motion are shown below. 
Throughout this thesis subscript 1 is for the towing (leading) ship and subscript 2 is for 


the towed (trailing) ship. 

(m 2 -X il2 )u 2 = -mv 2 r 2 -R 2 +T cos(^ 2 + y) (1) 

(m 2 - Y v2 )v 2 - Y. 2 r 2 = -m 2 r 2 u 2 + Y v2 v 2 + Y r2 r 2 - T sin (y/ 2 + y) (2) 

(1 zi ~ A^ 2)^2 ~-^v 2^2 ~N v2 v 2 + N r2 r 2 — Tx p2 sin(^ 2 +y) (3) 

(m, - Y n )Vj - Y n r x = -m x r x u x + Y vl v x + Y rX r x + T sin(^ + y) (4) 

(4i - ) ? 'i - = A vl v. + N rl r x - Tx pX sin(^ + y) (5) 

and 

Vi = r \ (6) 

y/ 2 = r 2 (7) 


The detailed explanation of the derivation of the equations of the motion, and the 
assumptions used can be found in the reference [1]. 

In these equations u is the surge velocity, v is sway velocity, r is the yaw velocity, 
R is the resistance of the vessel moving through body of water, and T is the tension in the 
connection (rope, cable, or some other mechanism) between the two towing ships. 

The kinematic relations are as shown below; 

y x =u j sin y/ x + v x cos y/ x (8) 

y 2 = u 2 sin y/ 2 + v 2 cos y/ 2 (9) 

If we define 
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y — yi ~y\ 


( 10 ) 


then 

y = y 2 -yi (U) 

These geometric relations are explained in Figure 3. 



Substituting equation (8) and equation (9) into equation (10), we get 
y = u 2 sin y/ 2 + v 2 cos y/ 2 -u l sin y x - v x cos y/ ] 

From the geometry of the figure we can see that; 
y 2 + x 2 sin - [y l - x , sin \j / x ) 

nn-i /IX - f r 


which can be rewritten as 


smy ■ 


v 1 ( 

- + sm \f /o +x 


/ / 


p i 



( 12 ) 


(13) 


(14) 
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B. LINEARIZATION 

In order to analyze the stability of our system, we should first linearize the 
system. During linearization we assumed that the velocities of the both vessels are 
constant and identical which in non-dimensional terms means 

u x = u 2 = 1 ( 15 ) 

Another assumption is that we are dealing with small heading angles, which gives 
us 

sin Wi,2=V'i,2 ( 16 ) 

and 

cos^j 2 = 1 (17) 

Since we are dealing with constant surge velocity, equation (1) is dropped, and we 
are left with the equations (2), (3), (4), and (5). 

When we substitute equations (15), (16), (17) into equations (2), (3), (4), (5), (12) 


and (13) we come up with 

(in 2 ~ Y n )v 2 - Y. 2 r 2 = -m 2 r 2 + Y v2 v 2 + Y r2 r 2 -T(y/ 2 +y) (18) 

(JZ2 — Nrl)n 2 — ^vl^l —N v 2 V 2 Y N r2 r 2 — Tx p2 (y/ 2 + y) ( 19 ) 

(/«! - Y. x )vj - Y n r x = -m x r x + K vl v, + Y rX r x + T(y/ X + y) (20) 

(4, - N n )fi - N n v x = N vl v x + N rX r x - Tx pX (if/, +y) (21) 

y = V2+v 2 ~W\- Vj (22) 

Y = ^+ l -[x p2 Y2+ x P x¥() ( 23 ) 

The above equations can be rewritten into a standard matrix form as follows: 

^2 — ^2vv V 2 ^2vr r 2 ^2v (^2 7 ) (24) 
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*2 i ^2rv’ ? 2 4,2+2 "*"^2r(f^2 + y) (25) 

'i = 4,'i + 5|„2- + «i, (y/, + f) (26) 

'1 = A„y, + '4.f + B, r (y/, + y) (27) 

After the necessary mathematical steps we can get the coefficients as follows 
4, = [(l, 2 -N tl )Y, 1 +(N, 1 Y n )h[(m 2 -Y a ll, 1 -N tI )-N il Y n \ (28) 

=[()', 2-»nl(V-,-V(J-(V(A;..)]e[(/«,-f,)(/.. N n Y n \ (29) 

4, =[-(2,2 -(V«)r-fc 2 Ji, 2 )J + [(m 2 -4X4 -N n )-N n Y n ] (30) 

=[(y, 1 n„ ! )+N jm, i,] (3i) 

4=fc-» 1 K 2 + 2V r2 (m 2 -rJMK -4X4 -AT f2 )-JV, 2 r, 2 ] (32) 

B 2r =\-TN ll -Tx p2 (m 1 -Y n h[(m 1 -Y n ll^-N n )-N n Y n \ (33) 

4, = [(4 -N,X, + (4,4 )M(™,~4 X/„ -N H )-N H T n ] (34) 

K, = [(4 -»>,X4 -4,)+(4,4)MK -4X4 -JV„)-JV„K„] (35) 

= 1(4 -A7,)7 - -te,J - iv)MK -4X4(36) 

4, =[(n.Af„)+Af„('», "4 )]+[(”>, -4X4-4,)-4,4] (37) 

Kr = 1(4 + 4, ('», - 4 )]+[('», - 4 X 4 -4,)-4,4l (38) 

s„ =[nv„ ] + 74 ,( 121 ,- 4 )]+[(in, - 4 X 4 -4)-44l (39) 

Now we can write down the complete system in a 7x7 matrix form as 
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^2 


X 2 



r 2 

A 


A 

fi 

= U] 

fi 

y 


y 

¥2 


¥2 



¥x_ 


where 


U] 


^2w ^2 w 0 0 

^2 rv A 2rr 0 0 

0 0 A lvv A lvr 

0 0 A lrv A hr 

10-10 
0 10 0 
0 0 0 1 


B 


2v 


B , 


B 2v + ' 


B-x 


2v p2 


B 2 r + 


1 

B 2r X p2 


/ 

B \r X p2 


l 

K K x 

i 

K 

i 

o 
o 
0 


P 2 


/ 


1 

0 

0 


B 2v X p\ 


B 2r X p\ 


K + 


B.x 


K+- 

-l 

o 

0 


hWpl 

/ 


/ 


(40) 


(41) 


Following linearization of the system, the next step is to find the eigenvalues of 
the A matrix. These will determine whether the system is stable or not. The matlab 
program shown in Appendix D was written and it was used to find the eigenvalues of the 
system. 

C. A ZERO EIGENVALUE 

Several different values for tension and for length of the connection between the 
two ships were tried using the program shown in Appendix D. In every instant, a zero 
eigenvalue was found. Therefore, it seemed impossible for the system to reach a stable 
condition. Simulations also showed the same effect. In order to test whether the existence 
of a zero eigenvalue was just a coincidence or an inherent property of the system, the 
matlab program shown in Appendix E is used. This performed a symbolic manipulation 
utilizing the Matlab symbolic manipulation toolbox. It was found that regardless of the 
numerical values of the geometric properties and the hydrodynamic coefficients of the 
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system, the characteristic equation of this system has always a zero constant, which 
proved the existence of a zero eigenvalue. Physically, this zero eigenvalue can be 
explained by the lack of any directional stabilizing effects on the leading ship. It can be 
easily verified, by observing the signs of the towline restoring forces and moments on the 
two ships that the towing force has a stabilizing effect on the trailing ship but a 
destabilizing effect on the leading ship. Since ships in the horizontal plane cannot exhibit 
directional stability, it is necessary in order to continue with the analysis to introduce 
some kind of directional control on the leading ship. This is analyzed in the following 
section. 

D. RUDDER CONTROL 

It was shown in the previous section that some kind of directional control on the 
leading ship is necessary to proceed with the towing analysis. We do this my means of a 
rudder control law. The equations of motion (4), (5) are changed because of inclusion of 


the rudder, but equations (2) and (3) are same. 

(*n 2 - Y v2 )v 2 - Yf 2 l% 2 = -m 2 r 2 u 2 + Y v2 v 2 + Y r2 r 2 - T sin (y/ 2 + y) (2) 

(/ Z2 — ^ i-2 )^2 ~^v2^2 ~ - kk v2 V 2 + ^ r 2^2 ~ p2 Y ) 0 ) 

(/«, - Y vt )v, - Y r[ r t = -m x r x u x + Y vX v x + Y rX r x + T sin(^ + y) + Y s S (42) 

(4i - N n~ N n v x = N vX v x + N rl r , - Tx pl sin(^ +y) + N s S (43) 

where 

8 = ~ k v V\ ~ Kyx - k A ~ k y y i ( 44 ) 


The rudder control law shown in equation (44) is a typical full state feedback 
control law. This was chosen as a representative control law and it by no means limits the 
applicability of the results that follow. The procedure would be identical even if a 
different rudder control law were employed. The feedback gains were calculated by using 
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the standard pole placement scripts in Matlab’s control system toolbox. In order to 
present the results in a compact form, we chose the control time constant as the 
representative selection criterion for the gains. Different choices do not alter the 
qualitative features of out results. The control time constant is defined as a non- 
dimensional response time for the system. It is measured in dimensionless seconds, with 
one dimensionless second being the time that it takes for a ship to travel one ship length. 
Typically, a system reaches steady state within three time constants. Thus, a high time 
constant signifies a rather slow leading ship control law, while a small time constant 
shows a fast (responsive) control. The closed loop poles are simply set as the negative 
inverse of the time constant. 

E. SAMPLE SIMULATION RESULTS 

We want to give one example to each stable and instable condition to clarify some 
of the concepts mentioned above. These results were obtained using the matlab code 
shown in Appendix H. 

1. Stable Condition 

The figures below are obtained for the following parameters; 

Tension=0.005 
Length=1.2 
Time Constant=l 

The results show a slow convergence to the nominal equilibrium state; i.e., 
straight line motion along the x-axis. These results will be confirmed with 
the stability analysis of the next chapter. 
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Figure 4.Ship Offset vs. Time 



Figure 5. Leading Ship Rudder Angle vs. Time 
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Figure 6. Towline Angle vs. Time 



Figure 7. Leading Ship Heading Angle vs. Time 
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Figure 8. Trailing Ship Heading Angle vs. Time 

2. Unstable Condition 

The figures shown below were obtained for the following parameters: 
Tension=0.005 
Length=1.2 
Time Constant=T .2 

The results clearly show an oscillatory divergence, and the system is 
unstable despite the stabilizing effect of the leading ship control law. The 
stability analysis of the following chapter will explain these results. 
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Figure 9. Ship Offsets vs. Time 



Figure 10. Leading Ship Rudder Angle vs. Time 
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Figure 11. Towline Angle vs. Time 



Figure 12. Leading Ship Fleading Angle vs. Time 
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Figure 13. Trailing Ship Heading Angle vs. Time 
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III. STABILITY ANALYSIS 


A. LINEARIZATION 

Since the equations of motion have changed by the inclusion effect of rudder, we 
have to again linearize the equations to see whether the system is stable or not. We will use 

the same assumptions that are mentioned before in Chapter 2 Section 2. 

When we linearize equations (2), (3), (8), (9), (13), (42), (43) we get 


( m 2 - Y v2 >2 - Y, 2 r 2 = -m 2 r 2 + Y v2 v 2 + Y r2 r 2 -T(y/ 2 +y) (45) 

{1 zi ~ Y ( " 2 )f 2 ~N i2 v 2 — Y v2 v 2 + Y ; . 2 r 2 — Tx p2 (l// 2 + y) (46) 

K - Y n )vj - Y n r x = -wy, + 7 vl Vj + Ky, + T(y/ X + y) + Y s S (47) 

( 4 i - # „ )r, - N n v, = N vl v, + N rA - Tx pl (ys l+ y) + N s S (48) 


y, = W\ + Vj 


(49) 


y =w2+ v 2 


(50) 


sin/ = 


y 2 +x p 2 y / 2 -U -x P iV\) 

i 


(51) 


We can convert the equations of motion into a matrix form in order to make it 
easier study the stability of the system. Since we changed only the equations of the first 
ship, equations (24), (25) are the same, but equations (26) and (27) are different 


^2 ^2 w r 2 ~^ B 2viV^2 7 ) (24) 

C — ^2rv V 2 -^2n- r 2 + B 2r^P2 7) (25) 

Vi = 4wh + 4 h y, + 5 lv + y) ■+ C h d (52) 

fi = 4rv V l + AnA + B lr iWl + ?) + C lr <? (53) 
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The coefficients are the same as in Chapter 2 Section 2, with the following 
exceptions: 


c„=[(/, -N„)Y S + AW.MK-n,X/ a -N„)-N*r„) (54) 

c ir = M (1 + N s (,n, - r„ )] -r [(,»,, - r H )(/„ - N„ ) - N„Y„ ] (55) 


As mentioned in Chapter 1, all hydrodynamic coefficients are taken from 
reference [2] with the exception of Y s and N s . The matlab program shown in Appendix 

F was used to calculate these coefficients. Since no data were available on the SLICE 
rudder design, we made an assumption of a ship turning radius of three ship lengths under 
fifteen degrees of rudder, and calculated the necessary values of the rudder hydrodynamic 
coefficients to achieve that turning radius. 

As mentioned before, the rudder feedback coefficients k ¥ ,k Y , k r and k y are 

calculated by standard pole placement techniques. The matlab program shown in 
Appendix G performs this calculation. 

Now we can set our matrix to analysis the stability of the system. The new system 
matrix is 8x8 and is written as shown below 


^2 


V 2 

c 


V 2 

A 


V 1 

K 

y 2 

= U] 

f 2 

Ti 


y i 

¥2 


¥2 

¥i\ 


¥x_ 


where 


(56) 
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0 


M= 


0 0 

0 0 

An - GA Aw — C\A, 
- C,k, 


A, 

o 

0 

0 

1 


-c,A 


A, 

/ 

i 

- 

A v X p2 

^2v X p\ 

l 

A, 

A r 

B 1 Bl ' Xpl 

Blr X p\ 

/ 

/ 


2r l 

l 

K 

-A 

- - C lv k y 

B\ v X p2 

B u . + KXp ' 

i 

l 


l 

n l 

K 

_A< 

--C k 

K*p2 

B + B '’ Xp ' 

l 

i 

( “ / l r K y 

l 

K+ l 

0 

0 


1 

0 

0 

0 


0 

1 

0 

0 


0 

0 

0 

0 


0 

0 


-CA 


~C lr k v 


(57) 


After fonning the matrix, we used the matlab code in Appendix H to find the 
eigenvalues of the matrix, which will dictate whether the system is stable or not. When 
we run the code, we experienced that the stability of the system changes according to 
different values of tension, length and the time constant. A systematic series of runs was 
conducted in order to ascertain the effects of these parameters. 


B. DEGREE OF STABILITY 

The critical eigenvalue of the 8x8 system; i.e., the eigenvalue that dictates 
stability or instability is the one that has the largest real part. We define this as the degree 
of stability of the system. A positive value indicates instability while a negative value 
indicates stability. A systematic series of runs was then conducted where the towline 
length, the towline tension, and the leading ship control time constant were varied. For 
this analysis we used the matlab code in Appendix I for varying length of the connection 
between the ships, the code in Appendix J for varying tension on the connection, and the 
code in Appendix K for varying control time constant. 

The results for different tow lengths are shown in Appendix A, for different 
tensions in Appendix B, and for different time constants in Appendix C. Based on these 
results we can draw the following conclusions: 

1. Relatively large time constants (a slow control law) cannot guarantee stability 
under towing. Sufficiently fast control laws may, depending on tension and 
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towline length, stabilize the complete system. Therefore, towing stability must be 
a consideration during leading ship control law design. 

2. Sufficiently large values of the tension (i.e.; trailing ship resistance) may stabilize 
the complete two-body system. This result is consistent with earlier observations 
reported in Ref. [3]. 

3. For very large values of the towline length, the effect of the leading ship control 
appears to be small. Therefore, the previously reported results in the literature [3] 
can be considered as a large towline length of the two-body model developed in 
this work. 


C. REGIONS OF STABILITY 

The previous results can be summarized in a single graph designating regions of 
stability and instability. This graph is shown in Figure 14 which was produced by 
utilizing the Matlab code shown in Appendix L. The graph shows the critical value of the 
control time constant for stability versus the towline length, and is parameterized by the 
tension in the towline. Combinations of (Tc,L,T) below the corresponding curve will 
yield stable response, while combinations that are located above the curve will result in 
system instability. A comparison between the values of the parameters used in the 
previous two simulation cases and the results shown in Figure 14, demonstrates the 
agreement between numerical integrations and the theoretical predictions of the behavior 
of the system based on our stability analysis. 

Finally, we present a comparison between our results and the results presented in 
Ref. [3], which were based on the 4 th order system. In reference [3], the first necessary 
stability criterion based on the trailing ship is 


x 


p 



(58) 


where x p is the non-dimensional length between the center of the ship and the 
towline connection point. 
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Different T Values 



Figure 14. The Stability Region 


The second necessary criterion is 


T>-°^~ (59) 

a x 

where T is the tension of the towed ship. If the above tension is negative then the 
criterion in equation (59) is inactive and is, therefore, automatically satisfied. 

The parameters in the above equations are: 


a 2 = ( B 0 C l D l -A 0 D l ~)+-(B q C 1 D 2 + B n C 2 D ] -2A 0 D { D 2 ) 
+— {b o C 2 D 2 — A 0 D 2 ) 

and 

a 2 = B 0 C 0 D ] +-(b 0 C 0 D 2 -B^eJ 


(60) 


(61) 


The coefficients used in equations (60), and (61) are as follows 
A 0 =(Y v -m)(N r -f z )-N v Y r (62) 
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(63) 


B 0 = (K - m)N r + Y v (N r -I z )+ N v (m -Y r )~ Y r N v 


Y v N r +N v {m-Y r ) 

(64) 

(m-Y.)x p + N v 

(65) 

-(N r -I z )+(m-Y v )x p 2 +x p {N v + Y r ) 

(66) 

N v ~Y v x p 

(67) 

-Y v x p 2 + N v x p + (Y r - 7, )x p - N r + TV, 

(68) 

0 

(69) 


Using the hydrodynamic coefficients of the Kaimalino estimated in [2], and 
presented in the Matlab code in several of the Appendices in this thesis, we can see that 
x is always greater than the ratio of N v to F . Furthermore, the critical tension is always 

negative as is estimated by the Matlab code shown in Appendix M. As a result, the 
previous stability criteria would have predicted a system that would always be stable and 
would have missed the stability and instability regions shown in Figure 14. Therefore, in 
towing stability studies, the designer should consider not only the effect of the trailing 
ship, but also the effect of the leading ship. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

This thesis presented a comprehensive study of the stability properties of the two- 
body ship towing problem. The results of this study can help in establishing rational 
guidelines which should be followed in order to ensure stability and safety of close 
proximity ship towing operations. The main results from this study can be summarized as 
follows: 

1. Relatively large time constants (a slow control law) cannot guarantee stability 
under towing. Sufficiently fast control laws may, depending on tension and 
towline length, stabilize the complete system. Therefore, towing stability must be 
a consideration during leading ship control law design. This is in contrast with 
previous studies where the dynamics of the leading ship were not considered 
important in the overall analysis. 

2. Sufficiently large values of the tension (i.e.; trailing ship resistance) may stabilize 
the complete two-body system. This result is consistent with earlier observations 
reported in the literature. 

3. For very large values of the towline length, the effect of the leading ship control 
becomes less pronounced. Therefore, previously reported results in the literature 
can be considered as a large towline length of the two-body model developed in 
this work. 
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B. RECOMMENDATIONS 

Recommendations for continuing studies are the following: 

1. Study the effect of additional geometrical parameters such as attachment point 
location on the regions of stability and instability. 

2. Perform a nonlinear analysis in order to analyze the mechanism of the loss of 
stability. Current results indicate that the predominant mechanism is through the 
generation of limit cycles (periodic solutions) but this needs to be substantiated 
further. 
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APPENDIX A. DEGREE OF STABILITY VS. TOWLINE LENGTH 


Constant T=0.001 



Fig A.l 

Constant T=0 003 



Fig A.2 
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Degree of Stability Degree of Stabilit 


Constant T=0.005 



Fig A.3 

Constant T=0.0075 



Fig A.4 
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Constant T=0.01 



Fig A.5 


Constant T c =0.1 



Fig A.6 
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Constant T c =0.3 



Fig A.7 


Constant T c =0 5 



Fig A.8 
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Constant T c =0.75 



Fig A.9 


Constant T c =1.0 



Fig A. 10 
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Constant T c =1.5 



Fig A. 11 


Constant T c =2.0 
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APPENDIX B. DEGREE OF STABILITY VS. TENSION 


Constant T c =0.1 



FigB.l 


Constant T c =0.3 



Fig B.2 
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Constant T c =0.5 



Fig B.3 


Constant T c =0.7 



Fig B.4 
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Constant T c =1.0 



Fig B.5 


Constant T c =1.25 



Fig B.6 
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Degree of Stability 


Constant T c =1.50 



Fig B.7 

Constant T c =1.75 



Fig B.8 
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Degree of Stability 


Constant T c =2.0 



Fig B.9 


Constant L=0.1 



Fig B. 10 
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Constant L=0.3 



Fig B. 11 


Constant L=0.5 



Fig B. 12 
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Constant L=0.7 



Fig B. 13 


Constant L=1.0 



Fig B.14 
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Constant L=1.25 



Fig B.15 


Constant L=1.50 



Fig B. 16 
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Constant L—1.75 



Fig B. 17 


Constant L=2.0 



Fig B.18 
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APPENDIX C. DEGREE OF STABILITY VS. TIME CONSTANT 


Constant T=0 001 



FigC.l 


Constant T=0.002 



FigC.2 
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Constant T=0.003 



FigC.3 


Constant T=0.004 



Fig C.4 
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Constant T=0.005 



Fig C.5 


Constant T=0 006 



Fig C.6 
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Fig C.7 


Constant T=0 008 



Fig C.8 
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Constant T=0.009 



Fig C.9 


Constant T=0.01 



Fig C. 10 
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Fig C. 11 
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Fig C. 13 
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Fig C. 14 
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Fig C.15 
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Fig C. 16 
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Fig C.18 
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Fig C. 19 
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APPENDIX D. MATLAB CODE FOR STABILITY ANALYSIS WITHOUT 

CONTROL 


%LTJG Mersin GOKCE 
% Thesis Work 

% Model for coupled stability analysis 
% Initialization 
ml=0.018078; m2=0.018; 
ul = l; 
u2 = l; 

T=l; 

L=l; 

xpl=0.5; xp2=0.5; 

lzl=0.0007; Iz2=0.00069412; 

Yvl=-0.07893; 

Yv2=-0.1183; 

Yrl=-0.004044; 

Yr2=-0.0042; 

Nvl=-0.016428; 

Nv2=-0.0187; 

Nrl=-0.010332; 

Nr2=-0.0176; 

Yvdotl=-0.051328; 

Yvdot2=-0.0184; 

Yrdot1=0.005617; 

Yrdot2=-0.0011; 

Nvdotl=-0.001945; 

Nvdot2=-0.0008489; 

Nrdotl=-0.00564; 

Nrdot2=-0.0090; 

A3=l/L; 

B3=xpl/L; 

C3=xp2/L; 

o, 

o 

A2vv=[((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

A2vr=[((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 
Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2v=[-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

A2rv=[(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

A2rr=[((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 
Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r=[-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

Alvv=[((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl-Yrdotl)]; 

Alvr=[(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blv=[((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alrv=[(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)] ; 

Alrr=[((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)]; 
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Blr=[(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)] ; 

A=zeros(7); 

A(1,1)=A2vv; 

A (1,2)=A2vr; 

A(1,5)=B2v/L; 

A(1, 6)=B2v+((B2v*xp2)/L); 

A(1,7)=(B2v*xpl)/L; 

A(2,1)=A2rv; 

A (2,2)=A2rr; 

A(2,5)=B2r/L; 

A(2,6)=B2r+((B2r*xp2)/L); 

A(2,7)=(B2r*xpl)/L; 

A(3, 3)=Alvv; 

A(3,4)=Alvr; 

A(3, 5)=Blv/L; 

A(3, 6) = (Blv*xp2)/L; 

A(3, 7)=Blv+( (Blv*xpl)/L); 

A(4,3)=Alrv; 

A(4,4)=Alrr; 

A (4,5)=Blr/L; 

A(4,6)=(Blr*xp2)/L; 

A(4,7)=Blr+ ( (Blr*xpl)/L); 

A(5, 1)=1; 

A(5, 3)=-l; 

A(5, 6)=1; 

A(5, 7)=-l; 

A(6, 2)=1; 

A(7,4)=1; 

eig(A) 


54 



APPENDIX E. MATLAB CODE TO PROVE ZERO EIGENVALUE 


syms al bl el fl gl 
syms a2 b2 e2 f2 g2 
syms c3 d3 e3 f3 g3 
syms c4 d4 e4 f4 g4 
syms h 

A=[al bl 0 0 el fl gl; 
a2 b2 0 0 e2 f2 g2; 

0 0 c3 d3 e3 f3 g3; 

0 0 c4 d4 e4 f4 g4; 

10-1001-1; 

0 1 0 0 0 0 0; 

0001000]; 

B=eye (7) .*h; 

C=det(A-B); 
h=0; 

D=subs(C); 

syms ml m2 ul u2 T L 

syms xpl xp2 Izl Iz2 

syms Yvl Yv2 Yrl Yr2 

syms Nvl Nv2 Nrl Nr2 

syms Yvdotl Yvdot2 Yrdotl Yrdot2 

syms Nvdotl Nvdot2 Nrdotl Nrdot2 

syms A2vv A2vr B2v 

syms A2rv A2rr B2r 

syms Alvv Alvr Blv 

syms Alrv Alrr Blr 

L=l; 

xpl=0.5; 
xp2=0.5; 

A2vv=[((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2-Nrdot2)) 
(Nvdot2*Yrdot2)] ; 

A2vr=[(Yr2-m2*u2)*(Iz2-Nrdot2)+(Nr2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 
Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2v=-[((Iz2-Nrdot2)*T)+(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2-Nrdot2)) 
(Nvdot2*Yrdot2)] ; 

A2rv=[(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

A2rr=[((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2 
Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r=-[(T*Nvdot2)+(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2))- 
(Nvdot2*Yrdot2)]; 

Alvv=[((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl)) 
(Nvdotl-Yrdotl)]; 

Alvr=[(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blv=[((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alrv=[(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alrr=[((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl 
Nrdotl))-(Nvdotl*Yrdotl)] ; 

Blr=[(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 


55 



al=A2vv; 
bl=A2vr; 
el=B2v/L; 

fl=B2v+((B2v*xp2)/L); 

gl=(B2v*xpl)/L; 

a2=A2rv; 

b2=A2rr; 

e2=B2r/L; 

f2=B2r+((B2r*xp2)/L); 

g2=(B2r*xpl)/L; 

c3=Alvv; 

d3=Alvr; 

e3=Blv/L; 

f3=(Blv*xp2)/L; 

g3=Blv+((Blv*xpl)/L); 

c4=Alrv; 

d4=Alrr; 

e4=Blr/L; 

f4=(Blr*xp2)/L; 

g4=Blr+((Blr*xpl)/L); 

E=subs(D) 
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APPENDIX F. MATLAB CODE TO FIND HYDRODYNAMIC COEFFICIENTS 


% to find Ydel and Ndel 
% Ndel=-0.5Ydel 
% R=1/(K*del) 

% K=[(Nvl*Ydel)-(Yvl*Ndel)]/[(Yvl*Nrl)-(Nvl*(Yrl-(ml*ul)))] 
% R=3 (Assumption,tactical diameter) 

% del=15 degree= 0.262 radian (Assumption,rudder) 
ml=.018078; 
ul=l ; 

Yvl=-0.07893; 

Yrl=-0.004044; 

Nvl=-0.016428; 

Nrl=-0.010332; 

R=3; 

del=-0.262 

Ydel=[Yvl*Nrl-(Nvl*(Yrl-ml))]/[(Nvl+(0.5*Yvl))*R*del] 
Ndel=-0.5*Ydel 
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APPENDIX G. MATLAB CODE TO FIND THE COEFFICIENTS OF DELTA 


% STEP 2 

% To Find Coefficients of Delta 
ml=.018078; 
ul = l; 

T=l; 

L=l; 

xpl=0.5; 

Izl=.0007; 

Yvl=-0.07893; 

Yrl=—0.004044; 

Nvl=-0.016428; 

Nrl=-0.010332; 

Yvdotl=-0.051328; 

Yrdot1=0.005617; 

Nvdotl=-0.001945; 

Nrdotl=-0.00564; 

Ydel=0.0103; 

Ndel=-0.0051; 


Alvv=[((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alvr=[(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)] ; 

Blv=[((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)] ; 

Clv=[((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)] ; 

Alrv=[(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Alrr=[((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 
Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blr=[(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)]; 

Clr=[(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl))- 
(Nvdotl*Yrdotl)] ; 

A=zeros(4); 

A(1,3)=1; 

A(2,2)=Alvv; 

A(2,3)=Alvr; 

A(3, 2)=Alrv; 

A(3, 3)=Alrr; 

A(4,1)=1; 

A(4,2)=1; 

B=zeros(4,1); 

B(2, 1) =Clv; 

B(3, 1)=Clr; 

poles=[-0.1 -0.11 -0.12 -0.13]; 
k=place(A,B,poles) 
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APPENDIX H. MATLAB CODE FOR STABILITY ANALYSIS AND 

SIMULATION WITH CONTROL 

% Main program 

% Eigenvalue analysis and numerical simulation 


% Parameters: 
% T 
% L 
% TC 
% xpl 
amidships) 

% xp2 

amidships) 

% DeltaT = 
% SimTime = 

o, 

o 


Nondimensional tension 

Towline length 

Control law time constant 

= Attachment point (leading ship - positive forward 

= Attachment point (trailing ship - positive aft 

Time step increment 
Simulation time 


T 

0 

.0075; 

L 

0 

.201; 

TC 

0 

• 6; 

bpole = 

- 

1/TC; 

xpl 

0 

• 5; 

xp2 = 

0 

• 5; 

DeltaT = 

0 

. 002; 

SimTime= 

2 

o 

o 

NPrint = 

1 

0; 


o, 

o 

o, 

o 

o, 

o 


Initial conditions for simulation 
yl,y2 must be consistent with L 


every 


(NPrint) 


points 


of 

of 


psil_old=0;vl_old=0;rl_old=0;yl_old=0.01; 
psi2_old=0;v2_old=0;r2_old=0;y2_old=0.00; 

o, 

o 


% Constants 


"6 


u 1 

= i; 

u2 

= 1; 

ml 

= 0.018078; 

m2 

= 0.018; 

Izl 

= 0.0007; 

Iz2 

= 0.00069412; 

Yvl 

= -0.07893; 

Yv2 

= -0.1183; 

Yrl 

= -0.004044; 

Yr2 

= -0.0042; 

Nvl 

= -0.016428; 

Nv2 

= -0.0187; 

Nr 1 

= -0.010332; 

Nr2 

= -0.0176; 

Yvdotl 

= -0.051328; 

Yvdot2 

= -0.0184; 

Yrdotl 

= 0.005617; 

Yrdot2 

= -0.0011; 

Nvdotl 

= -0.001945; 

Nvdot2 

= -0.0008489; 

Nrdotl 

= -0.00564; 

Nrdot2 

= -0.0090; 
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Ydel = 0.0103; 

Ndel = -0.0051; 

A3 = 1/L; 

B3 = xpl/L; 

C3 = xp2/L; 

D3 = -1/L; 

"o 

% Setup the matrix coefficients 

o, 

o 

A2vv = [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2-Nrdot2)) 

(Nvdot2*Yrdot2)] ; 

A2vr = [((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2-Yvdot2)*(Iz2 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2v = [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2-Nrdot2)) 

(Nvdot2*Yrdot2)] ; 

A2rv = [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2)) 

(Nvdot2*Yrdot2)]; 

A2rr = [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r = [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2-Nrdot2)) 

(Nvdot2*Yrdot2)] ; 

Alvv = [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl)) 

(Nvdotl*Yrdotl)] ; 

Alvr = [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl 

Nrdotl))-(Nvdotl*Yrdotl)] ; 

Blv = [((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl-Nrdotl)) 

(Nvdotl*Yrdotl)] ; 

Civ = [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)]/[((ml-Yvdotl)*(Izl-Nrdotl)) 

(Nvdotl*Yrdotl)] ; 

Alrv = [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl)) 

(Nvdotl*Yrdotl)]; 

Alrr = [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl 

Nrdotl))-(Nvdotl*Yrdotl)] ; 

Blr = [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/(((ml-Yvdotl)*(Izl-Nrdotl)) 

(Nvdotl*Yrdotl)] ; 

Clr = [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl-Nrdotl)) 

(Nvdotl*Yrdotl)] ; 

o, 

o 

% Find control gains 

o, 

o 


c 

= 

zeros (4); 

C (1,3) 

= 

i; 

c ( 2 , 2 ) 

= 

Alvv; 

C(2,3) 

= 

Alvr; 

C(3,2) 

= 

Alrv; 

C (3, 3) 

= 

Alrr; 

C (4,1) 

= 

i; 

C (4,2 ) 

= 

i; 

D 

= 

zeros (4,1); 

D(2,1) 

= 

Civ; 

D(3,1) 

= 

Clr; 

poles 

= 

[bpole bpole-0.05 bpole-0.10 bpole-0.15]; 

k 

= 

place (C,D,poles); 

Kpsi 

= 

k (1,1) ; 

Kv 

= 

k (1,2) ; 

Kr 

= 

k (1,3) ; 

Ky 

= 

k (1,4) ; 
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o, 

o 

% A matrix 
"6 

A = zeros (8); 

A (1,1) = A2vv; 

A (1,2) = A2vr; 

A (1,5) = B2v/L; 

A (1,6) = -B2v/L; 

A(l,7) = B2v+((B2v*xp2)/L); 

A (1,8) = (B2v*xpl)/L; 

A (2,1) = A2rv; 

A (2,2) = A2rr; 

A (2,5) = B2r/L; 

A (2,6) = -B2r/L; 

A(2,7) = B2r+((B2r*xp2)/L) ; 

A (2,8) = (B2r*xpl)/L; 

A (3,3) = Alvv-(Clv*Kv) ; 

A (3,4) = Alvr-(Clv*Kr) ; 

A(3,5) = Blv/L; 

A (3,6) = -Blv/L-(Clv*Ky); 

A(3,7) = (Blv*xp2)/L; 

A (3,8) = Blv+((Blv*xpl)/L)-(Clv*Kpsi); 

A (4,3) = Alrv-(Clr*Kv); 

A (4,4) = Alrr-(Clr*Kr); 

A (4,5) = Blr/L; 

A (4,6) = -Blr/L-(Clr*Ky); 

A (4,7) = (Blr*xp2)/L; 

A (4,8) = Blr+((Blr*xpl)/L)-(Clr*Kpsi); 

A (5, 1) = 1; 

A (5, 7) = 1; 

A (6, 3) = 1; 

A ( 6, 8 ) = 1; 

A (7,2 ) = 1; 

A (8,4 ) = 1; 

O, 

o 

% Find eigenvalues 
"6 

B=eig(A) 

o, 

o 

% Start simulation (Euler fixed step) 

o, 

o 

NT=SimTime/DeltaT; 
iPrint=l;iStore=l; 
for i=l:NT, 

o, 

o 

% Rudder angle 
"6 

delta = - ( Kpsi*psil_old + Kv*vl_old + Kr*rl_old + Ky*yl_old ); 

o, 

o 

% Limit rudder angle between -0.4 and +0.4 radians 

o, 

o 

if delta > 0.4 
delta=0.4; 

end 

if delta < -0.4 
delta=-0.4; 

end 
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% Calculate angle gamma 

o, 

o 

gamma=asin((y2_old+xp2*sin(psi2_old)-yl_old+xpl*sin(psil_old))/L); 

o, 

o 

% Calculate xdot=f(x) 

"5 

psildot = rl_old; 

vldot = Alvv*vl_old + Alvr*rl_old + Clv*delta + 

Blv*sin(gamma+psil_old); 

rldot = Alrr*rl_old + Alrv*vl_old + Clr*delta + 

Blr*sin(gamma+psil_old) ; 

yldot = ul*sin(psil_old) + vl_old*cos(psil_old); 
psi2dot = r2_old; 

v2dot = A2vv*v2_old + A2vr*r2_old + B2v*sin(gamma+psi2_old); 
r2dot = A2rr*r2_old + A2rv*v2_old + B2r*sin(gamma+psi2_old); 
y2dot = u2*sin(psi2_old) + v2_old*cos(psi2_old); 

o, 

o 

% Calculate new x = old x + Dt * xdot 


o, 

o 


psil_new 

= psil_old 

+ 

DeltaT 

•k 

psildot; 

vl_new 

= vl_old 

+ 

DeltaT 

-k 

vldot; 

rl_new 

= rl_old 

+ 

DeltaT 

•k 

rldot; 

yl_new 

= yl_old 

+ 

DeltaT 

•k 

yldot; 

psi2_new 

= psi2_old 

+ 

DeltaT 

•k 

psi2dot; 

v2_new 

= v2_old 

+ 

DeltaT 

•k 

v2dot; 

r2_new 

= r2_old 

+ 

DeltaT 

•k 

r2dot; 

y2_new 

= y2_old 

+ 

DeltaT 

•k 

y2dot; 


o, 

o 


% Ensure psi is between 2pi and -2pi 

o, 

o 

if psil_new > (2*pi) 

psil_new=psil_new-2*pi; 

end 

if psi2_new > (2*pi) 

psi2_new=psi2_new-2*pi; 

end 

if psil_new < (-2*pi) 

psil_new=psil_new+2*pi; 

end 

if psi2_new < (-2*pi) 

psil_new=psi2_new+2*pi; 

end 

o, 

o 

% Update x 

o, 

o 

psil_old=psil_new;vl_old=vl_new;rl_old=rl_new;yl_old=yl_new; 
psi2_old=psi2_new;v2_old=v2_new;r2_old=r2_new;y2_old=y2_new; 

o, 

o 

% Store results every NPrint time steps 

o, 

o 

iPrint=iPrint+l; 
if iPrint > NPrint 

iStore = iStore+1; 

time(iStore) = i*DeltaT; 
del (iStore) = delta; 
gam(iStore) = gamma; 
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psil (iStore) 

= psil_new; 

psi2 (iStore) 

= psi2_new; 

vl (iStore) 

= vl_new; 

v2(iStore) 

= v2_new; 

rl(iStore) 

= rl_new; 

r2(iStore) 

= r2_new; 

yl (iStore) 

= yl_new; 

y2(iStore) 

= y2_new; 

iPrint 

= i; 


end 

end 

"6 

% Results 

o, 

o 


figure (1) 

plot(time,yl,time,y2);legend('yl 
figure(2) 

plot(time,gam*57.2958);xlabel ( ' t 
figure(3) 

plot(time,de1*57.2958);xlabel('t 
figure (4) 

plot(time,psil*57.2958);xlabel ('t 
figure (5) 

plot(time,psi2*57.2958);xlabel ('t 


, 'y2');xlabel('y');ylabel('y');grid 
) ; ylabel('\gamma (deg) ');grid 


ylabel (' 

1 \delta 

(deg) ' ) 

; grid 

;ylabel 

( '\psi_l 

(deg)’ 

1 );grid 

;ylabel 

( '\psi_2 

(deg)’ 

1 );grid 
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APPENDIX I. MATLAB CODE TO FIND THE DEGREE OF STABILITY FOR 

DIFFERENT TOWLINE LENGTHS 

% Main program for Length 

% Eigenvalue analysis and numerical simulation 
"6 

% Parameters: 

% T = Nondimensional tension 

% L = Towline length 

% TC = Control law time constant 

% xpl = Attachment point (leading ship - positive forward of 

amidships) 

% xp2 = Attachment point (trailing ship - positive aft of 

amidships) 


% Constants 

o, 

o 


u 1 

= 

1; 

u2 

= 

i; 

ml 

= 

0.018078; 

m2 

= 

0.018; 

Izl 

= 

0.0007; 

Iz2 

= 

0.00069412; 

Yvl 

= 

-0.07893; 

Yv2 

= 

-0.1183; 

Yrl 

= 

-0.004044; 

Yr2 

= 

-0.0042; 

Nvl 

= 

-0.016428; 

Nv2 

= 

-0.0187; 

Nr 1 

= 

-0.010332; 

Nr2 

= 

-0.0176; 

Yvdotl 

= 

-0.051328; 

Yvdot2 

= 

-0.0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0.0011; 

Nvdotl 

= 

-0.001945; 

Nvdot2 

= 

-0.0008489; 

Nrdotl 

= 

-0.00564; 

Nrdot2 

= 

-0.0090; 

Ydel 

= 

0.0103; 

Ndel 

= 

-0.0051; 

T 

= 

\—1 
O 

o 

TC 

= 

2.0; 

bpole 

= 

-1/TC; 

xpl 

= 

0.5; 

xp2 

= 

0.5; 

index= 

0; 


for i 

=0 

.1:0.01:2.0; 

index 

=index+l; 

L 


= i; 

A3 


= 1/L; 

B3 


= xpl/L; 

C3 


= xp2/L; 

D3 


= -1/L; 


% Setup the matrix coefficients 

o, 

o 
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A2vv = [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2vr = [((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)]; 

B2v = [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rv = [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rr = [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r = [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

Alvv = [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl))-(Nvdotl*Yrdotl)] ; 

Alvr = [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blv = [((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)] ; 

Civ = [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrv = [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrr = [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml- 

Yvdotl )*(Izl-Nrdotl))-(Nvdotl*Yrdotl)]; 

Blr = [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Clr = [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

"6 

% Find control gains 

o, 

o 

C = zeros (4); 

C(1,3) = 1; 

C (2,2) = Alvv; 

C(2,3) = Alvr; 

C(3,2) = Alrv; 

C(3,3) = Alrr; 

C(4,1) = 1; 

C(4,2 ) = 1; 

D = zeros (4,1); 

D (2,1) = Civ; 

D(3,1) = Clr; 

poles = [bpole bpole-0.05 bpole-0.10 bpole-0.15]; 

k = place(C,D,poles); 

Kp s i = k(1,1) ; 

Kv = k (1,2 ) ; 

Kr = k (1,3) ; 

Ky = k (1,4) ; 

"6 

% A matrix 


A = zeros (8); 

A (1,1) = A2vv; 

A(1,2) = A2vr; 

A(1,5) = B2v/L; 

A (1,6) = -B2v/L; 

A(l,7) = B2v+( (B2v*xp2)/L); 
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A (1,8 ) 
A(2,1) 
A(2,2) 
A(2,5) 
A(2,6) 
A(2,7) 
A(2,8) 
A (3, 3) 
A (3, 4) 
A (3, 5) 
A (3, 6) 
A (3, 7) 
A (3, 8) 
A (4,3) 
A (4,4 ) 
A (4,5) 
A (4, 6) 
A (4,7 ) 
A (4,8 ) 
A (5, 1) 
A (5, 7) 
A (6, 3) 
A ( 6, 8 ) 
A (7,2 ) 
A (8,4 ) 

o, 

o 

% Find 

o, 

o 


= (B2v*xpl)/L; 

= A2rv; 

= A2rr; 

= B2r/L; 

= -B2r/L; 

= B2r+((B2r*xp2)/L); 

= (B2r*xpl)/L; 

= Alvv-(Clv*Kv); 

= Alvr-(Clv*Kr) ; 

= Blv/L; 

= -Blv/L-(Clv*Ky); 

= (Blv*xp2)/L; 

= Blv+((Blv*xpl)/L)-(Clv*Kpsi) ; 
= Alrv-(Clr*Kv) ; 

= Alrr-(Clr*Kr) ; 

= Blr/L; 

= -Blr/L-(Clr*Ky) ; 

= (Blr*xp2)/L; 

= Blr+((Blr*xpl)/L)-(Clr*Kpsi) ; 

= i; 

= 1 ; 

= 1 ; 

= i; 

= i; 

= i; 

eigenvalues 


B=eig(A); 

F=real(B); 

H(index)=max(F) ; 
Lv(index)=L; 

end 
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APPENDIX J. MATLAB CODE TO FIND THE DEGREE OF STABILITY FOR 

DIFFERENT TENSIONS 

% Main program for Tension 

% Eigenvalue analysis and numerical simulation 


% Parameters: 

% T = Nondimensional tension 

% L = Towline length 

% TC = Control law time constant 

% xpl = Attachment point (leading ship - positive forward of 

amidships) 

% xp2 = Attachment point (trailing ship - positive aft of 

amidships) 


% Constants 

o, 

o 


u 1 

= 

1; 


u2 

= 

i; 


ml 

= 

0 . 

018078; 

m2 

= 

0 . 

018; 

Izl 

= 

0 . 

0007; 

Iz2 

= 

o. 

00069412; 

Yvl 

= 

-o. 

07893; 

Yv2 

= 

-o. 

1183; 

Yrl 

= 

-o. 

004044; 

Yr2 

= 

-o. 

0042; 

Nvl 

= 

-o. 

016428; 

Nv2 

= 

-o. 

0187; 

Nr 1 

= 

-o. 

010332; 

Nr2 

= 

-o. 

0176; 

Yvdotl 

= 

-0 . 

051328; 

Yvdot2 

= 

-0 . 

0184; 

Yrdotl 

= 

0 . 

005617; 

Yrdot2 

= 

-0 . 

0011; 

Nvdotl 

= 

-0 . 

001945; 

Nvdot2 

= 

-0 . 

0008489; 

Nrdotl 

= 

-0 . 

00564; 

Nrdot2 

= 

-0 . 

0090; 

Ydel 

= 

0 . 

0103; 

Ndel 

= 

-0 . 

0051; 

L 

= 

2.0 

t 

TC 

= 

2.0 

r 

bpole 

= 

-1/ 

TC; 

xpl 

= 

0.5 

r 

xp2 

= 

0.5 

r 

A3 

= 

1/L; 

B3 

= 

xpl/L; 

C3 

= 

xp2/L; 

D3 

= 

-1/L; 

index= 

0; 



for i 

=0 

o 

o 

1:0.001:0 


index=index+l; 

T = i; 

% Setup the matrix coefficients 

o, 

o 
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A2vv = [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2vr = [((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)]; 

B2v = [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rv = [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rr = [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r = [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

Alvv = [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl))-(Nvdotl*Yrdotl)] ; 

Alvr = [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blv = [((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)] ; 

Civ = [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrv = [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrr = [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml- 

Yvdotl )*(Izl-Nrdotl))-(Nvdotl*Yrdotl)]; 

Blr = [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Clr = [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

"6 

% Find control gains 

o, 

o 

C = zeros (4); 

C(1,3) = 1; 

C (2,2) = Alvv; 

C(2,3) = Alvr; 

C(3,2) = Alrv; 

C(3,3) = Alrr; 

C(4,1) = 1; 

C(4,2 ) = 1; 

D = zeros (4,1); 

D (2,1) = Civ; 

D(3,1) = Clr; 

poles = [bpole bpole-0.05 bpole-0.10 bpole-0.15]; 

k = place(C,D,poles); 

Kp s i = k(1,1) ; 

Kv = k (1,2 ) ; 

Kr = k (1,3) ; 

Ky = k (1,4) ; 

"6 

% A matrix 


A = zeros (8); 

A (1,1) = A2vv; 

A(1,2) = A2vr; 

A(1,5) = B2v/L; 

A (1,6) = -B2v/L; 

A(l,7) = B2v+( (B2v*xp2)/L); 
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A (1,8 ) 
A(2,1) 
A(2,2) 
A(2,5) 
A(2,6) 
A(2,7) 
A(2,8) 
A (3, 3) 
A (3, 4) 
A (3, 5) 
A (3, 6) 
A (3, 7) 
A (3, 8) 
A (4,3) 
A (4,4 ) 
A (4,5) 
A (4, 6) 
A (4,7 ) 
A (4,8 ) 
A (5, 1) 
A (5, 7) 
A (6, 3) 
A ( 6, 8 ) 
A (7,2 ) 
A (8,4 ) 

o, 

o 

% Find 

o, 

o 


= (B2v*xpl)/L; 

= A2rv; 

= A2rr; 

= B2r/L; 

= -B2r/L; 

= B2r+((B2r*xp2)/L); 

= (B2r*xpl)/L; 

= Alvv-(Clv*Kv); 

= Alvr-(Clv*Kr) ; 

= Blv/L; 

= -Blv/L-(Clv*Ky); 

= (Blv*xp2)/L; 

= Blv+((Blv*xpl)/L)-(Clv*Kpsi) ; 
= Alrv-(Clr*Kv) ; 

= Alrr-(Clr*Kr) ; 

= Blr/L; 

= -Blr/L-(Clr*Ky) ; 

= (Blr*xp2)/L; 

= Blr+((Blr*xpl)/L)-(Clr*Kpsi) ; 

= i; 

= 1 ; 

= 1 ; 

= i; 

= i; 

= i; 

eigenvalues 


B=eig(A); 

F=real(B); 

H(index)=max(F) ; 
Tv(index)=T; 

end 
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APPENDIX K. MATLAB CODE TO FIND THE DEGREE OF STABILITY FOR 

DIFFERENT TIME CONSTANTS 

% Main program Time Constant 

% Eigenvalue analysis and numerical simulation 


% Parameters: 
% T 
% L 
% TC 
% xpl 
amidships) 

% xp2 

amidships) 


Nondimensional tension 

Towline length 

Control law time constant 

= Attachment point (leading ship - positive forward 
= Attachment point (trailing ship - positive aft 


of 

of 


% Constants 

o, 

o 


u 1 

= 

1 

r 

u2 

= 

1 

r 

ml 

= 

0 

. 018078; 

m2 

= 

0 

.018; 

Izl 

= 

0 

.0007; 

Iz2 

= 

0 

.00069412; 

Yvl 

= 

-0 

.07893; 

Yv2 

= 

-0 

.1183; 

Yrl 

= 

-0 

.004044; 

Yr2 

= 

-0 

.0042; 

Nvl 

= 

-0 

.016428; 

Nv2 

= 

-0 

.0187; 

Nr 1 

= 

-0 

.010332; 

Nr2 

= 

-0 

.0176; 

Yvdotl 

= 

-0 

.051328; 

Yvdot2 

= 

-0 

.0184; 

Yrdotl 

= 

0 

. 005617; 

Yrdot2 

= 

-0 

.0011; 

Nvdotl 

= 

-0 

. 001945; 

Nvdot2 

= 

-0 

. 0008489; 

Nrdotl 

= 

-0 

.00564; 

Nrdot2 

= 

-0 

.0090; 

Ydel 

= 

0 

.0103; 

Ndel 

= 

-0 

.0051; 

T 

= 

0 . 

\—1 
o 

L 

= 

2 . 

0; 

xpl 

= 

0 . 

5; 

xp2 

= 

0 . 

5; 

A3 

= 

1 

/L; 

B3 

= 

xpl/L; 

C3 

= 

xp2/L; 

D3 

= 

-1 

/ L; 

index= 

0; 



for i 

=0 

. 1 

:0.01:2.0; 


index=index+l; 

TC =i; 

bpole = -1/TC; 

% Setup the matrix coefficients 

o, 

o 
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A2vv = [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2vr = [((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)]; 

B2v = [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rv = [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rr = [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r = [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

Alvv = [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl))-(Nvdotl*Yrdotl)] ; 

Alvr = [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Blv = [((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)] ; 

Civ = [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrv = [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrr = [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml- 

Yvdotl )*(Izl-Nrdotl))-(Nvdotl*Yrdotl)]; 

Blr = [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Clr = [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

"6 

% Find control gains 

o, 

o 

C = zeros (4); 

C(1,3) = 1; 

C (2,2) = Alvv; 

C(2,3) = Alvr; 

C(3,2) = Alrv; 

C(3,3) = Alrr; 

C(4,1) = 1; 

C(4,2 ) = 1; 

D = zeros (4,1); 

D (2,1) = Civ; 

D(3,1) = Clr; 

poles = [bpole bpole-0.05 bpole-0.10 bpole-0.15]; 

k = place(C,D,poles); 

Kp s i = k(1,1) ; 

Kv = k (1,2 ) ; 

Kr = k (1,3) ; 

Ky = k (1,4) ; 

"6 

% A matrix 


A = zeros (8); 

A (1,1) = A2vv; 

A(1,2) = A2vr; 

A(1,5) = B2v/L; 

A (1,6) = -B2v/L; 

A(l,7) = B2v+( (B2v*xp2)/L); 
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A (1,8 ) 
A(2,1) 
A(2,2) 
A(2,5) 
A(2,6) 
A(2,7) 
A(2,8) 
A (3, 3) 
A (3, 4) 
A (3, 5) 
A (3, 6) 
A (3, 7) 
A (3, 8) 
A (4,3) 
A (4,4 ) 
A (4,5) 
A (4, 6) 
A (4,7 ) 
A (4,8 ) 
A (5, 1) 
A (5, 7) 
A (6, 3) 
A ( 6, 8 ) 
A (7,2 ) 
A (8,4 ) 

o, 

o 

% Find 

o, 

o 


= (B2v*xpl)/L; 

= A2rv; 

= A2rr; 

= B2r/L; 

= -B2r/L; 

= B2r+((B2r*xp2)/L); 

= (B2r*xpl)/L; 

= Alvv-(Clv*Kv); 

= Alvr-(Clv*Kr) ; 

= Blv/L; 

= -Blv/L-(Clv*Ky); 

= (Blv*xp2)/L; 

= Blv+((Blv*xpl)/L)-(Clv*Kpsi) ; 
= Alrv-(Clr*Kv) ; 

= Alrr-(Clr*Kr) ; 

= Blr/L; 

= -Blr/L-(Clr*Ky) ; 

= (Blr*xp2)/L; 

= Blr+((Blr*xpl)/L)-(Clr*Kpsi) ; 

= i; 

= 1 ; 

= 1 ; 

= i; 

= i; 

= i; 

eigenvalues 


B=eig(A); 

F=real(B); 

H(index)=max(F) ; 
TCv(index)=TC; 

end 
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APPENDIX L. MATLAB CODE TO FIND THE REGION OF STABILITY 


% Stability graph - TC vs. L - constant T 

o_ 

o 

% Parameters: 


% T 
% L 
% TC 
% xpl 
amidships) 
% xp2 

amidships) 


= Nondimensional tension 
= Towline length 
= Control law time constant 

= Attachment point (leading ship - positive forward 

= Attachment point (trailing ship - positive aft 


of 

of 


% Constants 

o, 

o 


u 1 

= 

i; 

u2 

= 

i; 

ml 

= 

0.018078; 

m2 

= 

0.018; 

Izl 

= 

0.0007; 

Iz2 

= 

0.00069412; 

Yvl 

= 

-0.07893; 

Yv2 

= 

-0.1183; 

Yrl 

= 

-0.004044; 

Yr2 

= 

-0.0042; 

Nvl 

= 

-0.016428; 

Nv2 

= 

-0.0187; 

Nr 1 

= 

-0.010332; 

Nr2 

= 

-0.0176; 

Yvdotl 

= 

-0.051328; 

Yvdot2 

= 

-0.0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0.0011; 

Nvdotl 

= 

-0.001945; 

Nvdot2 

= 

-0.0008489; 

Nrdotl 

= 

-0.00564; 

Nrdot2 

= 

-0.0090; 

Ydel 

= 

0.0103; 

Ndel 

= 

-0.0051; 

indexL 

= 

0; 

index 

= 

0; 

xpl 

= 

0.5; 

xp2 

= 

0.5; 


o_ 

o 


% Enter constant T 


T = input('T = ' ); 

o, 

o 

% Start loop on length 

"o 

for iL =0.1:0.01:2.0; 

indexL = indexL+1; 

indexTC = 0; 

L = iL 


L_v(indexL)= L; 

A3 = 1/L; 
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B3 = xpl/L; 

C3 = xp2/L; 

D3 = -1/L; 

o, 

o 

% Loop on TC 

o, 

o 

for iTC =0.1:0.01:2.0; 
indexTC=indexTC+l; 

TC =iTC; 

bpole = -1/TC; 

TC_v(indexTC)=TC; 

"6 

% Setup the matrix coefficients 

o, 

o 

A2vv = [((Iz2-Nrdot2)*Yv2)+(Nv2*Yrdot2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2vr = [((Yr2-(m2*u2))*(Iz2-Nrdot2))+(Nr2*Yrdot2)]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)]; 

B2v = [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2~ 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rv = [(Yv2*Nvdot2)+(Nv2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

A2rr = [((Yr2-(m2*u2))*Nvdot2)+(Nr2*(m2-Yvdot2))]/[((m2- 

Yvdot2)*(Iz2-Nrdot2))-(Nvdot2*Yrdot2)] ; 

B2r = [-(T*Nvdot2)-(T*xp2*(m2-Yvdot2))]/[((m2-Yvdot2)*(Iz2- 

Nrdot2))-(Nvdot2*Yrdot2)] ; 

Alvv = [((Izl-Nrdotl)*Yvl)+(Nvl*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl))-(Nvdotl*Yrdotl)] ; 

Alvr = [(Yrl-ml*ul)*(Izl-Nrdotl)+(Nrl*Yrdotl)]/[((ml- 

Yvdotl) *(Izl-Nrdotl))-(Nvdotl*Yrdotl)] ; 

Blv = [((Izl-Nrdotl)*T)-(Yrdotl*T*xpl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)] ; 

Civ = [((Izl-Nrdotl)*Ydel)+(Ndel*Yrdotl)]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

Alrv = [(Yvl*Nvdotl)+(Nvl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) ) -(Nvdotl*Yrdotl)]; 

Alrr = [((Yrl-(ml*ul))*Nvdotl)+(Nrl*(ml-Yvdotl))]/[((ml- 

Yvdotl )*(Izl-Nrdotl))-(Nvdotl*Yrdotl)]; 

Blr = [(T*Nvdotl)-(T*xpl*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)] ; 

Clr = [(Ydel*Nvdotl)+(Ndel*(ml-Yvdotl))]/[((ml-Yvdotl)*(Izl- 

Nrdotl) )-(Nvdotl*Yrdotl)]; 

o, 

o 

% Find control gains 

o, 

o 


c 

= zeros(4); 

C (1,3) 

= i; 

C(2,2) 

= Alvv; 

C(2,3) 

= Alvr; 

C(3,2) 

= Alrv; 

C(3,3) 

= Alrr; 

C (4,1) 

= i; 

C (4,2) 

= 1; 

D 

= zeros(4,1); 

D(2,1) 

= Civ; 

D(3,1) 

= Clr; 

poles 

= [bpole bpole-0.05 bpole-0.10 bpole-0.15]; 
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k = place(C,D,poles); 

Kp s i = k(1, 1) ; 

Kv = k (1,2 ) ; 

Kr = k (1,3) ; 

Ky = k (1,4) ; 

"6 

% A matrix 
"6 

A = zeros (8); 

A (1,1) = A2vv; 

A(1,2) = A2vr; 

A(1,5) = B2v/L; 

A (1,6) = -B2v/L; 

A(l,7) = B2v+((B2v*xp2)/L); 

A (1,8) = (B2v*xpl)/L; 

A (2,1) = A2rv; 

A (2,2) = A2rr; 

A (2,5) = B2r/L; 

A (2,6) = -B2r/L; 

A(2,7) = B2r+((B2r*xp2)/L) ; 

A (2,8) = (B2r*xpl)/L; 

A (3,3) = Alvv-(Clv*Kv) ; 

A (3,4) = Alvr-(Clv*Kr) ; 

A (3,5) = Blv/L; 

A (3,6) = -Blv/L-(Clv*Ky); 

A(3,7) = (Blv*xp2)/L; 

A (3,8) = Blv+((Blv*xpl)/L)-(Clv*Kpsi); 

A (4,3) = Alrv-(Clr*Kv) ; 

A (4,4) = Alrr-(Clr*Kr) ; 

A (4,5) = Blr/L; 

A (4,6) = -Blr/L-(Clr*Ky) ; 

A (4,7) = (Blr*xp2)/L; 

A (4,8) = Blr+ ( (Blr*xpl)/L)-(Clr*Kpsi); 

A (5, 1) = 1; 

A (5, 7) = 1; 

A (6, 3) = 1; 

A(6,8 ) = 1; 

A (7,2) = 1; 

A (8,4 ) = 1; 

O, 

o 

% Find eigenvalues 
"6 

B=eig(A); 

F=real(B); 

H_new=max (F); 

o, 

o 

% Detect if H changes its sign 
"6 

if indexTC==l 

H_old=H_new; 

else 

PROD=H_new*H_old; 
if PROD > 0 
"6 

% H does not change sign - keep going 
"6 

H_old=H_new; 
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else 


% H changed its sign - find critical TC by 

interpolation 

"6 

index=index+l; 

TC_crit(index)=-(H_old*TC_v(indexTC)- 
H_new*TC_v(indexTC-1))/(H_new-H_old); 

H_old=H_new; 

end 

end 

end 

end 

o, 

o 

% Plot (L, crit TC) curve 

o, 

o 

%plot(L_v, TC_crit) , xlabel('L'),ylabel('{T_C}_{\rm 
crit}'),title ( 'Constant T'),grid 


linear 
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APPENDIX M. MATLAB CODE TO EVALUATE THE STABILITY CRITERIA 

BASED ON ONLY TRAILING SHIP 

%Check the Stability Criteria in the paper 

%For the Second Ship 

%First Criteria xp2 > Nv2/Yv2 

%Second Criteria T > -(Alfa2/Alfal) 

m2=.018; 

u2 = l; 

T=0.01; 

L=1.8541; 
xp2=0.5; 

Iz2=.00069412; 

Yv2=-0.1183; 

Yr2=-0.0042; 

Nv2=-0.0187; 

Nr2=-0.0176; 

Yvdot2=-0.0184; 

Yrdot2=-0.0011; 

Nvdot2=-0.0008489; 

Nrdot2=-0.0090; 

Xp2=(Nv2/Yv2) 

A=[((Yvdot2-m2)*(Nrdot2-Iz2))-(Nvdot2*Yrdot2)]; %A0 
B=[((Yvdot2-m2)*Nr2)+(Yv2*(Nrdot2-Iz2))+(Nvdot2*(m2-Yr2))- 
(Yrdot2*Nv2)]; %B0 

C=[(Yv2*Nr2)+(Nv2*(m2-Yr2))]; %C0 
Cl=[ ((m2-Yvdot2)*xp2)+Nvdot2] ; 

C2=[-(Nrdot2-Iz2) + ((m2-Yvdot2)*(xp2 A 2)) + (xp2*(Nvdot2+Yrdot2)) ] ; 

Dl=[Nv2-(Yv2*xp2)]; 

D2=[-(Yv2*(xp2 A 2)) + (Nv2*xp2) + ((Yr2-Yvdot2)*xp2)-Nr2+Nvdot2] ; 

E1=0; 

Alfal=[ ( (B*C1*D1)-(A*(D1 A 2) ) ) + ( (1/L)* ( (B*C1*D2) + (B*C2*D1)~ 

(2*A*D1*D2) ) ) + ( (1/ (L A 2) )* ( (B*C2*D2)-(A*(D2 A 2) ) ) ) ] 

Alfa2=[ (B*C*D1) + ((1/L)*((B*C*D2)- (B A 2*E1))) ] 

Tl=—(Alfa2/Alfal) 
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